#include "outline/grand_poly_bezier_builder2.h"
#include "outline/bezier.h"
#include "outline/bezier_guide.h"
#include <limits>

namespace outline {

length_record length_table::operator [] (uint l) const
{
	//position of last node of open spline
	//uint last_min = size + 1 - max_estimated_length;
	uint last_max = size - 1;

	length_record r;
	r.min_start = 0;//std::max<int>(0, last_min - l + 1);
	r.max_start = last_max - l + 1;
	assert (r.min_start <= r.max_start);
	r.can_be_closed = l + max_estimated_length - 1 > size;
	return r;
}

void junction_table::reset(uint size, uint max_estimated_length)
{
	vec_.clear();
	min_jcount_ = ~0;
	if (size < 4)
		return;
	for (uint j = 0;; ++j) {
		junction_record r;
		r.min_len = 3*j + 4;
		r.max_len = std::min((j + 1)*(max_estimated_length - 1) + 1, size - 2);
		if (r.min_len > r.max_len)
			break;
		r.can_be_closed = r.max_len + max_estimated_length - 1 > size;
		if (r.can_be_closed && min_jcount_ == ~0)
			min_jcount_ = j;
		vec_.push_back(r);
	}
}

template <typename T>
inline void bottom_limit(T & v, T limit)
	{ if (v < limit) v = limit; }

template <typename T>
inline void top_limit(T & v, T limit)
	{ if (v > limit) v = limit; }

template <typename T>
inline void limit(T & v, T bottom, T top)
	{ bottom_limit(v, bottom); top_limit(v, top); }

template <typename R>
void grand_params<R>::put_within_the_limits()
{
	bottom_limit<R>(sharpness_penalty, 0);
	limit(low_max_iterations, 1U, 100000U);
	limit<R>(low_stop_mean_discrepancy, 0, 10);
	limit<R>(low_oversample_factor, 0.01, 10);
	limit(hi_max_iterations, 1U, 100000U);
	limit<R>(hi_stop_mean_discrepancy, 0, 10);
	limit<R>(hi_oversample_factor, 0.1, 10);
	limit(max_junctions, 1U, 100U);
	limit(min_junctions, 1U, max_junctions);
	bottom_limit<R>(penalty_for_each_junction, 0);
}

template <typename T, typename R>
void grand_poly_bezier_builder2<T, R>::reset(nodes_type const & polyline)
{
	reset_nodes(polyline);
	reset_estimations();
}

template <typename T, typename R>
void grand_poly_bezier_builder2<T, R>::reset_estimations()
{
	max_estimated_length_ = 3;
	max_estimated_start_ = 0;
	estimations_.clear();
	estimations_.resize(size_);

	length_table_.reset(size_, max_estimated_length_);
	junction_table_.reset(size_, max_estimated_length_);

	candidate_.clear();

	reset_assemblings();
}

template <typename T, typename R>
void grand_poly_bezier_builder2<T, R>::reset_assemblings()
{
	jvec_.clear();
}

template <typename T, typename R>
void grand_poly_bezier_builder2<T, R>::set_params(grand_params<R> nw)
{
	nw.put_within_the_limits();
	if (params_.low != nw.low || nw.low && (
			params_.low_max_iterations != nw.low_max_iterations ||
			params_.low_stop_mean_discrepancy != nw.low_stop_mean_discrepancy ||
			params_.low_oversample_factor != nw.low_oversample_factor
	))
		reset_estimations();
	else if (params_.penalty_for_each_junction != nw.penalty_for_each_junction ||
			 params_.sharpness_penalty != nw.sharpness_penalty)
		reset_assemblings();

	params_ = nw;
}

template <typename T, typename R>
void grand_poly_bezier_builder2<T, R>::estimate_more()
{
	if (all_lengths_estimated())
		return;

	reset_assemblings();

	++max_estimated_length_;
	bezier_type b;
	for (uint i = 0; i < size_; ++i) {
		R e = bezier_guide(
				nodes_, i, i + max_estimated_length_,
				params_.low_max_iterations,
				params_.low_stop_mean_discrepancy,
				params_.low_oversample_factor,
				false, b
			);
		assert (e >= 0);
		estimations_[i].push_back(estimation_type(e, 
			b.start_control() - b.start(), b.finish_control() - b.finish()));
		assert(estimations_[i].size() == max_estimated_length_ - 3);
	}

	length_table_.reset(size_, max_estimated_length_);
	junction_table_.reset(size_, max_estimated_length_);
}

/*template <typename T, typename R>
void grand_poly_bezier_builder2<T, R>::estimate_all()
{
	if (all_lengths_estimated())
		return;

	reset_estimations();

	for (uint i = 0; i < size_; ++i) {
		estimate_one(i);
	}

	max_estimated_length_ = size_ + 1;
	length_table_.reset(size_, max_estimated_length_);
	junction_table_.reset(size_, max_estimated_length_);
}*/

template <typename T, typename R>
void grand_poly_bezier_builder2<T, R>::estimate_one()
{
	if (all_starts_estimated())
		return;

	uint i = max_estimated_start_;
	estimations_[i].reserve(size_+1 - 3);

	incremental_bezier_guide<T,R> ibg;
	ibg.max_iterations = params_.low_max_iterations;
	ibg.stop_mean_discrepancy = params_.low_stop_mean_discrepancy;
	ibg.oversample_factor = params_.low_oversample_factor;

	//initialization: piece of minimal length
	R e = ibg.initial_fit(nodes_, i, i + 4);
	assert (e >= 0);
	estimations_[i].push_back(estimation_type(e, ibg.b.start_d(), ibg.b.finish_d()));
	assert(estimations_[i].size() == 1);

	//iterative incrementation
	for (uint l = 5; l <= size_ + 1; ++l) {
		R e = ibg.incremental_fit(i + l);
		assert (e >= 0);
		estimations_[i].push_back(estimation_type(e, ibg.b.start_d(), ibg.b.finish_d()));
		assert(estimations_[i].size() == l - 3);
	}

	++max_estimated_start_;
	if (all_starts_estimated()) {
		max_estimated_length_ = size_ + 1;
		length_table_.reset(size_, max_estimated_length_);
		junction_table_.reset(size_, max_estimated_length_);

		/*uint min_piece_len = 10;
		std::vector<partition> partitions(min_piece_len-1);
		uint begin;
		uint best_begin;
		R best_partition_d = std::numeric_limits<R>::max();
		for (begin = 0; begin < min_piece_len-1; ++begin) {
			closed_partition(begin, min_piece_len, partitions[begin]);
			R d = partitions[begin][0].d;
			assert (d >= 0);
			if (d < best_partition_d) {
				best_partition_d = d;
				best_begin = begin;
			}
		}
		assert(best_partition_d != std::numeric_limits<R>::max());
		partition const & p = partitions[best_begin];

		for (uint i = 0; i < p.size(); i = p[i].next) {
			candidate_.push_back(best_begin + i);
		}*/
	}
}

template <typename T, typename R>
void grand_poly_bezier_builder2<T,R>::assemble_more()
{
	if (all_junctions_assembled())
		return;

	if (jvec_.empty()) {
		init_jvec();
		return;
	}

	uint j = jvec_.size();
	unijunction_routes & uj = append_jvec();
	//unijunction_routes const & prev_uj = *jvec_[j-1];
	junction_record const & jr = junction_table_[j];
	for (uint l = jr.min_len; l <= jr.max_len; ++l) 
	{
		unilength_routes & unilen = uj[l - jr.min_len];
		length_record lr = length_table_[l];
		unilen.resize(lr.scount());
		for (uint s = lr.min_start; s <= lr.max_start; ++s) 
		{
			route & r = unilen[s - lr.min_start];
			r.d = std::numeric_limits<R>::max();
			uint const max_body_len = junction_table_[j-1].max_len;
			uint const min_body_len = junction_table_[j-1].min_len;
			uint min_head_len = std::max<int>(4, l - max_body_len + 1);
			assert (l > min_body_len);
			uint max_head_len = std::min(max_estimated_length_, l - min_body_len + 1);
			for (uint l0 = min_head_len; l0 <= max_head_len; ++l0) 
			{
				estimation_type const & head = estimations_[s][l0 - 4];
				assert (head.d >= 0);
				uint body_start = s + l0 - 1;
				uint body_len = l - l0 + 1;
				route const & body = get_jvec_at(j-1, body_len, body_start);
				assert (body.d >= 0);
				R d = head.d + body.d;
				if (d >= r.d)
					continue;
				uint neck_end = body.second;
				estimation_type const & neck = estimations_[body_start][neck_end - body_start + 1 - 4];
				if (params_.sharpness_penalty > 0)
					d += sharpness_penalty(params_.sharpness_penalty, head.d2, neck.d1);
				if (d < r.d) {
					r.d = d;
					r.second = body_start;
				}
			}
			assert (r.d < std::numeric_limits<R>::max());
		}
	}
}

template <typename T, typename R>
void grand_poly_bezier_builder2<T,R>::init_jvec()
{
	unijunction_routes & j0 = append_jvec();
	junction_record const & jr = junction_table_[0];
	for (uint l = jr.min_len; l <= jr.max_len; ++l) {
		unilength_routes & unilen = j0[l - jr.min_len];
		length_record lr = length_table_[l];
		unilen.resize(lr.scount());
		for (uint s = lr.min_start; s <= lr.max_start; ++s) {
			route & r = unilen[s - lr.min_start];
			r.d = estimations_[s][l - 4].d;
			assert (r.d >= 0);
			r.second = s + l - 1;
		}
	}
}

template <typename T, typename R>
typename grand_poly_bezier_builder2<T,R>::unijunction_routes & 
grand_poly_bezier_builder2<T,R>::append_jvec()
{
	uint j = jvec_.size();
	jvec_.push_back(unijunction_routes_sptr(new unijunction_routes(junction_table_[j].lcount())));
	return *jvec_.back();
}

template <typename T, typename R>
R grand_poly_bezier_builder2<T,R>::select_best(index_vector_type & ivec)
{
	if (!params_.low) {
		uint mj = std::min(params_.max_junctions, size_/3);
		ivec.resize(mj + 1);
		for (uint j = 0; j < mj; ++j) //uniformly distributed junctions
			ivec[j] = size_ * j / mj;
		ivec[mj] = ivec[0] + size_;
		return 0;
	}

	ivec.clear();
	if (jvec_.empty() || jvec_.size() < junction_table_.min_jcount()+1)
		return -1;

	assert (junction_table_[jvec_.size() - 1].can_be_closed);
	R best_d = std::numeric_limits<R>::max();
	uint best_j, best_l, best_b; ///< best junctions count, length and beginning of open spline

	assert (params_.min_junctions >= 1);
	if (max_estimated_length_ == size_ + 1 && (params_.min_junctions == 1 || jvec_.empty())) {
		for (uint b = 0; b < size_; ++b) {
			estimation_type const & e = estimations_[b][max_estimated_length_-4];
			R d = e.d;
			assert (d >= 0);
			if (params_.sharpness_penalty > 0)
				d += sharpness_penalty(params_.sharpness_penalty, e.d2, e.d1);
			if (d < best_d) {
				best_d = d;
				best_b = b;
				best_l = size_ + 1;
			}
		}
	}

	uint min_j = 0;
	if (params_.min_junctions > 2)
		min_j = std::min(params_.min_junctions-2, jvec_.size()-1);
	for (uint j = min_j; j < jvec_.size(); ++j) 
	{
		junction_record const & jr = junction_table_[j];
		if (!jr.can_be_closed)
			continue;
		R penalty = (j+1) * params_.penalty_for_each_junction;
		if (best_d < std::numeric_limits<R>::max()) {
			if (best_d < penalty)
				break;
			if (j+1 >= params_.max_junctions) //it's limit for the junctions in closed spline
				break;
		}
		unijunction_routes const & unijun = *jvec_[j];

		for (uint l = jr.min_len; l <= jr.max_len; ++l) 
		{
			length_record lr = length_table_[l];
			if (!lr.can_be_closed)
				continue;
			unilength_routes const & unilen = unijun[l - jr.min_len];

			for (uint b = lr.min_start; b <= lr.max_start; ++b) 
			{
				route const & r = unilen[b - lr.min_start];
				R d = r.d;
				assert(d >= 0);
				estimation_type const & closure = inline_get_estimation(b + l - 1, size_ + 2 - l);
				assert (closure.d >= 0);
				d += closure.d;
				d += penalty;
				if (d >= best_d)
					continue;
				if (params_.sharpness_penalty > 0) {
					uint lhead = r.second - b + 1;
					estimation_type const & head = estimations_[b][lhead - 4];
					estimation_type const * tail;
					if (j == 0)
						tail = &head;
					else {
						//find route's tail
						uint j0 = j-1;
						uint b0 = r.second;
						uint l0 = l - lhead + 1;
						while (j0 > 0) {
							route const & sub = get_jvec_at(j0, l0, b0);
							--j0;
							l0 -= sub.second - b0;
							b0 = sub.second;
							assert (b + l == b0 + l0);
						}
						tail = &estimations_[b0][l0 - 4];
					}
					d += sharpness_penalty(
						params_.sharpness_penalty, closure.d2, head.d1);
					d += sharpness_penalty(
						params_.sharpness_penalty, tail->d2, closure.d1);
				}
				if (d < best_d) {
					best_d = d;
					best_j = j;
					best_l = l;
					best_b = b;
				}
			}
		}
	}

	assert (best_d < std::numeric_limits<R>::max());

	//write down the best route
	ivec.clear();

	if (best_l == size_ + 1) {
		ivec.resize(2);
		ivec[0] = best_b;
		ivec[1] = best_b + size_;
		return best_d;
	}

	ivec.reserve(best_j + 3);
	ivec.push_back(best_b);

	uint j0 = best_j;
	uint b0 = best_b;
	uint l0 = best_l;
	while (j0 > 0) {
		route const & r = get_jvec_at(j0, l0, b0);
		--j0;
		l0 -= r.second - b0;
		b0 = r.second;
		ivec.push_back(b0);
		assert (best_b + best_l == b0 + l0);
	}

	ivec.push_back(b0 + l0 - 1);
	ivec.push_back(best_b + size_);

	return best_d;
}

template <typename T, typename R>
typename grand_poly_bezier_builder2<T,R>::estimation_type const & 
grand_poly_bezier_builder2<T,R>::get_estimation(uint begin, uint len, uint) const
{
	return inline_get_estimation(begin, len);
}

template <typename T, typename R>
void grand_poly_bezier_builder2<T,R>::get_junction_types(
	index_vector_type const & ivec, std::vector<junc_type> & junc) const
{
	if (!params_.low) {
		junc_type t = (params_.sharpness_penalty > 0) ? junc_smooth : junc_cusp;
		for (uint i = 0; i != junc.size(); ++i)
			junc[i] = t;
		return;
	}

	uint jcount = junc.size();
	assert(ivec.size() == jcount + 1);

	estimation_type const * e1 = &inline_get_estimation(ivec[jcount-1] % size_, ivec[jcount] - ivec[jcount-1] + 1);

	for (uint i = 0; i < jcount; ++i) {
		estimation_type const * e2 = &inline_get_estimation(ivec[i] % size_, ivec[i+1] - ivec[i] + 1);

		//if (square(e1->d2 + e2->d1) < 0.1 * (square(e1->d2) + square(e2->d1)))
		//	junc[i] = junc_symmetric;
		//else if (sharpness_penalty(real_type(1), e1->d2, e2->d1) < real_type(0.25)) //is angle between tangents less than 30 degrees
		//	junc[i] = junc_smooth;

		if (sharpness_penalty(real_type(1), e1->d2, e2->d1) < real_type(0.25)) //is angle between tangents less than 30 degrees
			junc[i] = junc_smooth;
		else
			junc[i] = junc_cusp;

		e1 = e2;
	}
}

template <typename T, typename R>
void grand_poly_bezier_builder2<T,R>::closed_partition(uint begin, uint min_len, partition & p)
{
	assert(all_starts_estimated());
	assert(min_len >= 4);
	uint max_len = 2*(min_len-1);

	p.resize(size_ - min_len + 2);
	uint i;
	for (i = min_len; i <= max_len; ++i) {
		p[size_-i+1] = step(estimations_[(begin + size_ - i + 1) % size_][i-4].d, ~0);
		assert (p[size_-i+1].d >= 0);
	}

	i = size_ - max_len;
	for (;;) {
		assert(begin + i < size_);
		p[i].d = std::numeric_limits<R>::max();
		for (uint l = min_len; l <= max_len; ++l) {
			uint next = i + l - 1;
			if (next >= p.size())
				break;
			assert (estimations_[begin + i][l-4].d >= 0);
			assert (p[next].d >= 0);
			R d = estimations_[begin + i][l-4].d + p[next].d;
			if (d < p[i].d)
				p[i] = step(d, next);
		}
		assert(p[i].d != std::numeric_limits<R>::max());
		if (i == 0)
			break;
		--i;
	}
}

template grand_poly_bezier_builder2<unsigned short>;
template grand_poly_bezier_builder2<unsigned short, double>;
template grand_poly_bezier_builder2<short>;
template grand_poly_bezier_builder2<int>;
template grand_poly_bezier_builder2<float>;
template grand_poly_bezier_builder2<double>;

} //namespace outline
